Measuring the elements of the optical density matrix 
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Most methods for experimentally reconstructing the quantum state of light involve determining 
a quasiprobability distribution such as the Wigner function. In this paper we present a scheme 
for measuring individual density matrix elements in the photon number state representation. Re- 
markably, the scheme is simple, involving two beam splitters and a reference field in a coherent 
state. 
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I. INTRODUCTION 

It is now well established that the quantum state of 
light can be measured. The first experimental evidence 
of this |l ] followed the work of Vogel and Risken , where 
it was shown that the Wigner function could be recon- 
structed from a complete ensemble of measured quadra- 
ture amplitude distributions. The authors of p] mea- 
sured the quadrature distributions using balanced homo- 
dyne techniques. In the case of inefficient homodyne de- 
tectors, a more general s-parameterized quasiprobability 
distribution is obtained resulting in a smoothed Wigner 
function. In either case, to obtain the quasiprobability 
phase space distribution from the measured data a rather 
complicated inverse transformation is required. 

Novel techniques which avoid this transformation are 
aimed at measuring the quasiprobability distribution 
more directly. This can be achieved, for example, in un- 
balanced homodyne counting experiments |3[ where 
a weighted sum of photocount statistics are combined to 
obtain a single point in the phase space distribution. The 
entire distribution is then obtained by scanning the mag- 
nitude and phase of the local oscillator over the region 
of interest while repeating the photon counting at each 
point. Perhaps the most direct method of obtaining a 
quasiprobability distribution is to use heterodyne or 
double homodyne Q |7j detection techniques where the 
Q function is measured. The Q function is related to the 
Wigner function through a convolution with a gaussian 
distribution which effectively washes out many of the in- 
teresting quantum features. It is possible to recover these 
features by deconvoluting the Q function, however this 
requires multiplying by an exponentially increasing func- 
tion thereby introducing a crucial dependence on sam- 
pling noise |^ . Other non-tomographic state reconstruc- 
tion schemes are proposed in P| for fields in a cavity, 
in for trapped atoms and in jl^ where use of a 
Schrodinger-cat state probe is suggested. Further dis- 
cussion of such techniques can be found in the recent 
review by Welsch and Vogel [|2| , the book by Leonhardt 



|l| and in Q. 

A different approach has been suggested by Steuer- 
nagel and Vaccaro , who have proposed an interesting 
nonrecursive scheme to measure not the quasiprobability 
distribution, but rather the density operator in the pho- 
ton number basis. The scheme is relatively direct in that 
only a finite number of different measurements are re- 
quired to determine each matrix element. However, the 
major disadvantage of this scheme is that determination 
of each matrix element, pmn, requires the preparation 
of a two-state superposition of |A^) and \M) for use as a 
probe field. Not only has the preparation of such fields 
not yet been achieved, but also the experiment requires 
a change of the probe field for the measurement of each 
matrix element. 

In this paper we investigate an alternative non- 
recursive scheme to that of Steuernagel and Vaccaro that 
also measures the individual density matrix elements in 
the photon number basis. Remarkably we find that this 
can be achieved with a single reference field that can be 
in an easily prepared coherent state. The only changes 
to the reference field needed to measure all the matrix 
elements are simple phase shifts. The integral transfor- 
mation required in the tomographic technique is avoided 
and is replaced by a summation of only four easily mea- 
surable probabilities. We find that this technique is par- 
ticulary suited for measurements of low intensity states, 
such as that used in where it offers some simplifica- 
tions over the tomographic methods. 



II. MEASUREMENT TECHNIQUE 

The density matrix element pmn in the number state 
representation of a density operator p is given by 



PMN^Tl{p\N){M\). 



(2.1) 



This can be compared with the probability for an out- 
come event e of a measurement on a system in state p, 
which is given from general quantum measurement the- 
ory by 



P(e)=Tr[/5n(e)], 



(2.2) 
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where n(e) is the element of a probability operator mea- 
sure (POM) associated with the event e. Comparing 
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FIG. 1: Apparatus for measuring the density matrix elements 
of light. BSl and BS2 are beam splitters. The field to be 
measured and a reference field in a coherent state are in the 
input modes c and a of BSl and BS2 respectively. A vacuum 
is in the input mode b and photon counters are in the output 
modes. Phase shifter PS adjusts the phase of the coherent 
state. 



and (2.2) suggests that if we could find a POM ele- 



ment equal to the operator |7V) (M| then we could find the 
matrix element pmn simply by measuring the probabil- 
ity P{e). This of course is not possible as the probability 
must be between zero and one but the matrix element 
need not even be real. However if we could synthesize 
the operator |iV)(A/| by a linear combination of different 
POM elements then we could find the matrix element 
from the same linear combination of the associated mea- 
surable probabilities. We adopt this operator synthesis 
approach in this paper. 

The proposed measurement technique uses the ar- 
rangement shown in Fig. |]. It consists of two symmetric 
beam splitters labelled BSl with input modes b and c 
and BS2 with input modes b and a. BS2 is a 50/50 beam 
splitter but we keep the transmission to reflection coeffi- 
cient ratio for BSl general for now. In the output modes 
are photodetectors Da, Dh and Dc- The input fields in 
modes b and a are respectively in a vacuum state |0)b and 
a coherent state \a)a- The optimum value of the ampli- 
tude of the coherent state will be discussed later. The 
field to be measured, in state pc, is in the input mode 
c of BSl. In the entry port of BS2 is a phase shifter 
PS capable of altering the phase of the coherent state, 
thereby changing the argument of the coherent state am- 
plitude. We let the amplitude of the coherent state be 
a = \a\ expiiip) at the entry of BS2, that is the argument 
(fi oi a incorporates the phase shift. We let (p he a func- 
tion of two numbers /3 and j, that is ip — ip{f3,j), that 
will be specified later. 

For simplicity, we assume that the distance between 
the beam splitters is an integer number of wavelengths 
of the light, which allows us to ignore the evolution of the 
light, which is just a phase shift, between the beam split- 
ters. If this is difficult experimentally, the discrepancy 
can be offset by an adjustment of the phase shifter. The 
complete time evolution operator is then R2R1 where i?2 



and i?i are the unitary operators for the action of beam 
splitters BS2 and BSl. Let e — {na,nh,nc) be the event 
that photodetectors Da, Db and Dc register ria, ni, and ric 
photocounts respectively. The probability for this event 



IS 



P/5j(e) = Trabc{R2Rl^ RlRl\na)aa{na\ <^ \nb)bb{nb\ 

®|7^c)ccK|), (2.3) 

where we have written the combined density operator for 
the three input fields as 



\a)aa{a\ (g) \0)bb{0\(g> Pc 



(2.4) 



and, as the subscripts imply, the trace is over the state 
spaces for all t hree modes. The subscript /3j on the prob- 
ability in (2.3) is to show explicitly that the probability 
is a function of the argument Lp[[3,j) of a, that is, it is 
a function of the setting of the phase shifter. Using the 
cyclic property of the trace we can write (p. 31) as 



where 



Pp,{e) = TYc /5cn^j(e) , (2.5) 



\lp,{e) - \q)cc{q\ (2.6) 



\q)c = b{m\ a{a\W,\na)a\nb)b\nc)c. (2.7) 



with 



From ( p. 51) we see that n^j(e) is an element of the POM 
for the measuring device that comprises all of the ar- 
rangement depicted in Fig. |l| except for the field to be 
measured. In general the elements of T\pj{e) are not nec- 
essarily orthogonal in that Ilpj{e)Iipj{e') is not necessar- 
ily zero for e 7^ e'. The origin of the non-orthogonality 
is the introduction of the two reference modes a and b. 
The effect of these two ancillary modes is to cube the 
dimensionality of the system space. In considering the 
measuring apparatus to consist of everything in Fig. 1 
except the state to be measured, we effectively reduce 
the apparatus to a single mode measuring device with 
many more POM elements than the dimensionality of 
the single mode. This means that the POM elements 
cannot all be orthogonal to one another. 

Our aim is to find a linear combination of POM ele- 
ments equal to the operator |iV)(M|. It is convenient to 
write this operator as \N) (TV -t- A| and consider separately 
the cases where A is even and odd. We examine first the 
case where A is even. Consider the particular detection 
event ei = (A/2, A/2, TV) in which the photodetectors Da, 
Db and Dc detect A/2, A/2 and N photocounts respec- 
tively. As shown in the Appendix, this turns out to be 
the optimum detection event for A even. The unitary 
operator Ri for BSl is given by fl^ 



i?i = exp[i?7(c^6 + ^^c)]. 



(2.8) 
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where cos rj = t and sin rj = r are the transmission and re- 
flection coefficients of BSl. i?2 for the 50/50 beam spUt- 
ter BS2 is a similar function of a and h with — 7r/4. 
Using these expressions with n^, and nj, equal to A/2, 
A/2 and N , we find from the Appendix that (2.7) be- 
comes 

\q)c= fneM-i{N + \-n)ip{P,j)]\n), (2.9) 



where 

In 



(2^)-V2^JV|c^|^+^-»(-^r)"-^^/^ 



exp(|a|V2)[(A + N - n)/2\\[{n - iV)/2]!VM 

(2.10) 

if n - is even and /„ = if n - iV is odd. The POM 
element (p.6[) for the detection event ei then becomes 



n^i(ei) 



N+X 

E 

n,m—N 



/n/m exp [?(n - m)(^(/3, j)]|n)cc(m|. 

(2.11) 



The terms with n — m odd are all zero. 



For A = we find from ( |2.1lD that the POM element in 
( p.6D is just proportional to \N)cc{N\, allowing us to find 
the diagonal elements pnn of the density matrix from the 
probability of detecting the event (0, 0, N). 

For A even and non-zero, we let (p{P,j) take particular 
values 



^(/3,J) = - + — 



(2.12) 



and consider a modified measurement procedure in which 
f3 is held constant but the value of j is cycled so that it 
takes all the integer values from to (A/2) — 1 with equal 
probability. This measurement procedure will have its 
own probability operator measure comprised of elements 
(e) . This will be different from our previous POM with 
elements n^j(e) because it describes a different measure- 
ment process. The POM element for detecting the event 
ei by means of this cycling procedure will be given by 

A/2-1 



n/3(ei) = y n«(ei) 



= J J2 {/„/;; exp[z(n-m)(/37r/A)] 

n,m—N 
A/2-1 

X ^ exp [z(ri — m)27rj/A]|n)cc(™|}- (2.13) 

The associated probability can be obtained in practice 
from the occurrence frequency of the event ei as we cy- 
cle through the values of j with the experiment being 
repeated an equal number of times for each value of j. 
Because we need only consider terms in (2.13) for which 



n — m is even, we can take the factor involving the sum- 
mation over j as zero unless n — m is zero or ±A, in which 
case it equals A/ 2. This gives us 

n/3(ei) = E \fn\'\n)cc{n\ + [fNrN+x 

X exp(-i7r/3)|Ar)cc(A^ + A| + h.c]{2.U) 

By choosing different values for /3, we obtain different 
cycling experiments, each with its own POM. Experimen- 
tally this means cycling through a diffe rent set of phase 
settings. It is not difficult to see from ( ^.14 ) that a lin- 
ear combination of POM elements 11^ (ei) with (3 taking 
the values 0, 1, 1/2 and 3/2 is required to synthesize the 
operator I Af)(A^ + A|. Specifically, 



|7V)(Ar + A| = 

no(ei) - ni(ei) + I ni/2(ei) - n3/2(ei) 



(2.15) 



N+X 



where 



1/2 



JnTn+x - t^''{^r|2f\aoal\ [^^^^ ^ (2.16) 

is a normalisation constant with o„ — {n\a). 

Taking the trace of the product of the density oper- 
ator pc of the field to be measured with both sides of 



15) gives the desired matrix element in terms of the 
measurable probabilities Pp{ei) for detecting the event 
(ei): 



{N + \\p,\N) = 

Po(ei) - Pi(ei) + z [Pi/2(ei) - ^3/2(61)] 



(2.17) 



JV+A 



for A non-zero and even. The complex conjugate of ( 2.17 ) 
is {N\p,\N + \). 

To find the density matrix element for A odd, we con- 
sider the detection event 62 = [(A + l)/2, (A - l)/2, A], 
which is shown in the Appendix to be the optimum de- 
tection event for this case. A derivation similar to that 
above eventually yields the associated POM element of 
the form 



N+X 

n/3j(e2) = Y 9n9m^'^v[i{n-m)ip{l3,j)]\n)cc{m\. 

71,771 — N 

(2.18) 

In ( 2.1^ ) n — m takes all integer values from +A to —A. 
We consider a measurement procedure where (3 is held 
constant but j takes all values from to A — 1 with equal 
probability. The POM element Iif}{e2) for detecting the 
event 62 with this procedure will be given by 



A-l 



n/5(e2) 



1 

3=0 



n«(e2) 
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^ X J2 {5«5mexp[i(n-TO)(/37r/A)] 

n,m—N 
A-1 

X ^ exp [i{n — m)2Trj / X]\n)cc{m\}. (2.19) 

The factor involving the summation over j is zero unless 
n — m = or ±A and is then equal to A. We find that the 
formula for the density matrix element in terms of the 
measurable probabilities Pp{e2) for detecting the event 
62 is, for A odd, 

{N + X\p,\N) = 

Po(e2) - Pi (62) + I [Pi/2(e2) - P3/2(e2)] ^2 20) 



where 



9N9*N+x = it'^'^{ir/2)^\aQal\ 



(A-l)/2 



N ) 

(2.21) 

We note the number of phase settings required for each 
experiment for the odd-A case is about twice that re- 
quired for the even- A case with a similar value of A. 
The number of experiments needed for the odd-A case 
can, however, be reduced by a factor of two as follows. 
While measuring the probability for the event 62, we 
can also measure the probability for another event 63 
= [(A - l)/2, (A + l)/2,iV]. It is possible to show that 
the probability ^^3(63) is equal to Pf3+i{e2)- Thus we 
only need two experiments with different value s of (3 to 
obtain all four terms in the numerator of (2.20). 



III. SOME PRACTICAL CONSIDERATIONS 

To illustrate the viability of our proposal we shall, in 
this section, show the extent to which physical imperfec- 
tions, such as a noisy local oscillator, detector inefficiency 
and dark counts, can be ignored or compensated for in a 
practical experiment. 

A more general measurement scheme would have an 
arbitrarily mixed reference state at the input of mode a 
in Fig. 1. That is \a)a, with a = |a|exp(i(^) at the 
entry of BS2, would be replaced by a density opera- 
tor exp(iAfa(^)/3a exp(— liVaV?)- This could represent, as 
a specific example, a noisy local oscillator. Following the 
derivation outlined above, it is straightforward to show 
the eff ect o f this is to replace the term |aoa^| in equa- 
tions ( ^.16 ) and ( ^.21 ) with the density matrix element 
(0 1 Pa I A). This is interesting because it shows how the 
density matrix element {N+X\pc\N) of the unknown field 
can be obtained directly from the density matrix element 
(0|/5q|A) of the reference field. Thus, provided we know 
the noise characteristics of the reference state, we do not 
require it to be a noiseless coherent state, or indeed any 
particular pure state, to use it to find the density matrix 



TABLE I: Truncated density matrix for a coherent state with 
a mean photon number of 0.5. 



0.6065 


0.4289 


0.2145 


0.0870 


0.0336 


0.4289 


0.3033 


0.1517 


0.0615 


0.0238 


0.2145 


0.1517 


0.0759 


0.0308 


0.0119 


0.0870 


0.0615 


0.0308 


0.0125 


0.0048 


0.0336 


0.0238 


0.0119 


0.0048 


0.0019 


TABLE IL Simulation of measured density matrix for a co- 
herent state with a mean photon number of 0.5 and detector 
inefficiency ri — 0.9 


0.6592 


0.4195 


0.1888 


0.0692 


0.0220 


0.4195 


0.2967 


0.1335 


0.0489 


0.0161 


0.1888 


0.1335 


0.0668 


0.0244 


0.0081 


0.0692 


0.0489 


0.0244 


0.0100 


0.0033 


0.0220 


0.0161 


0.0081 


0.0033 


0.0013 



elements of the unknown field. So we find in general that 
a noisy local oscillator can easily be used in the mea- 
surement scheme. A problem arises, however, if (0|/5a|A) 
is vanishingiy small in that the measured probabilities 
will c oincide with rare events as indicated by ( 2.17D and 
( ^.20D . This is the case when phase diffusion in the lo- 
cal oscillator is prominent, effectively diagonalizing the 
density matrix pa and removing all phase information. 
This can be avoided if both the reference field and the 
measured field, pc, are derived from a common source, 
a technique commonly exploited in experiments of this 
kind. 

Another practical issue concerns the extent to which 
inefficient photodetectors degrade the reliability of the 
measured data. The effect of inefficient photodetectors 
is to make the measurement process uncertain. For a 
given detector efficiency rj, the probability of detecting n 
photons, Pniv)^ is related to the probability of detecting 
m photons with a perfect detector, Pm(l)j by [Eoj 



Pniv) 



E 

m=Q 



^"(l-r/)"p„+™(l). (3.1) 



To illustrate what effect inefficient photodetectors have 
on the outcome of the experiment, some numerical calcu- 
lations were performed for a low intensity coherent input 
state with a mean photon number of 0.5. An example of 
the results are summarized in Tables | and 0, where the 
measured density matrix is displayed for a detector effi- 
ciency of 77 = 0.9. The reference state used in each simu- 
lation was a coherent state with a mean photon number 
of |ap = 0.5. As expected, as the efficiency decreases the 
relative error in the individual matrix elements increases. 
Fortunately it is possible to invert equation (^^) through 
a Bernoulli transformation and recover the exact proba- 
bilities from the detection statistics with sufficiently good 
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detectors 19'. This would allow accurate reconstruction 
of the density matrix. 

In addition, for weak fields in the quantum regime, 
with sufficiently long gating times, dead times need not 
be significant. If dead times are significant, more sophis- 
ticated detection methods are required for photon num- 
ber discrimination, such as replacing each detector with 
a multiport device such as described in [po[ . 

So far we have not specified the value of \a\ or t/r. 
The optimum values of these should maximize the de- 
nominators of ( |2.17 ) and ( 2.20| ), thereby avoiding quo- 
tients of small numbers. We find that the optimum value 
of |ap is A/2 and that of {t/r)"^ is 2N/X. As these are 
optimum values only, they need not be changed for the 
measurement of each matrix element and a reasonable 
compromise value should suffice, for example, for weak 
fields where the spread of values of N and A is not large. 

While the method proposed in this paper can be used 
to measure any individual density matrix element, it is 
not necessary to perform the same number of cycling ex- 
periments as matrix elements to find the density matrix. 
The matrix elements {N+X\pc\N) and their complex con- 
jugates for all values of N can be found from the same 
four cycling experiments. Also many phase settings can 
be used as parts of different cycling experiments, allow- 
ing further efficiencies. For example the setting ip{(3,j) 
— tt/2 can be used for /3 = 0, j = 1, A = 4 and (3 = 0, 
j = 2, A = 8 as well as for /? = 1, j = 0, A = 2 and so on. 



IV. CONCLUSION 

In this paper, we have extended the method of projec- 
tion synthesis in which a projector is synthesized by 
use of an exotic reference state, to a more general tech- 
nique of operator synthesis in which an operator is syn- 
thesized by a linear combination of POM elements. This 
provides a nonrecursive method for measuring individual 
density matrix elements of a light field. Remarkably, the 
technique is reasonably simple, involving only two beam 
splitters and a reference field which can be in an easily- 
prepared coherent state. In particular, for states that can 
be represented in a finite dimensional Hilbert space, this 
technique appears simpler than the tomographic methods 
in that only a finite number of different measurements 
are required to ascertain the complete density matrix. 
We have shown how detector inefficiency can be allowed 
for and have considered the effect of noise in the local os- 
cillator. We found that the local oscillator noise can be 
readily accounted for provided we know the correspond- 



ing mixed state description of the local oscillator. Inter- 
estingly, our method allows the density matrix elements 
of the unknown field to be obtained quite simply from the 
density matrix elements of a noisy local oscillator field, 
even when the unknown field is in a pure state. 
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APPENDIX 

In this appendix we derive the general form of \q)c 
defined by 



\q)c = b{0\Rl a{a\Rl\na)a\'nb)b\nc)c 



(A.l) 



With the unitary operator of a beam-splitter given by 



(A.2) 



where a and h are the annihilation operators for the in- 
put field modes, it can be shown that the beam-splitter 
transforms the corresponding creation operators, and 
, and the double mode vacuum according to 



R) a)R = ta^ ~irb^ 

i?t|o)fc|o), = |o)fc|o)e. 



(A.3) 
(A.4) 
(A.5) 



where t and r are the transmission and refiection co- 
efficients of the beam-splitter. In the case of BS2, a 
50/50 beam-splitter, t = r = 1/V2. By writing \na)a 
as (a^""/VnJ)|0)a, and similarly for we obtain 



Rl\na)a\nb)b 



2(na+nb)/2^„^!^j,! 



|0)a|0)b (A.6) 



and thus 



a{a\W2\na)a\nb)b\N)c 

{a* -ib^Y-{b^ -ia*)"^" 
~ 2("»+"'>)/2 exp(|a|2 /2)^/^^{Jm 



:|0),|iV)c (A.7) 



Writin g \N)^&s (ct^/VM)|0)c and using an equivalent 
form of (A.3-A.5), we obtain 



R\a{a\Rl\na)a\nb)b\N), 



i{tb^ 



'{tc 



rb^ 



N 



2(«„+n,)/2 exp(|a|2/2)\/«a!«6!iV! 



■|0).|0). 



(A., 
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where we have left the transmission and reflection coef- 
ficients of BSl as t and r. FinaUy, projecting onto the 
vacuum state in mode b gives us 



(-i)"H^[a* - rct]"» [a* + rct]"^ 
2(«„+«,)/2 exp(|a|2/2)V7IJ7I^ 

m=N 



where X — Ha + ni,. The specific notation 



\N)c 



(A.9) 



frn = gm(A/2, A/2) exp[i(iV + A - m)v3] 
ffm = <Zm[(A + l)/2, (A - l)/2] exp[i(iV + A - m)#.10) 

is used in the text. The exphcit form of qm{na, ^b) is not 
actuaUy needed. What is important is an expression for 



qN{na,nb)q*i<[j^\{nainh)- This can be derived from ( |A.9| ) 
by evaluating qNina^ni,) and gjv+A("'a: "■b) separately to 
give 

QNina, nb)qlf+xi'i^a, rib) 
= {-ir^aoaie^ir/2)^{l){^+'f\ (A.U) 

where a„ = {n\a). For a mixed reference state with den- 
sity operator exp{i Naip)pa exp(— zA^gy?) at the entry of 
BS2, aoa\ in (A. 11) is replaced by (0|/5o|A) exp(— lAi/?). 



It is not difficult to sec that the modulus of (A.ll) is 
maximized when Ua = A/2 if A is even and w hen Ug — 
(A± l)/2 if A is odd. Thus the quotients in (2.17) and 
( ^■20|) will have optimum numerators and denominators 
for the values of Ua and ni, we have used in this paper. 
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